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We study the two-dimensional bond-diluted XY and six-state clock models by Monte Carlo simu- 
lation with cluster spin updates. Various concentrations of depleted bonds were simulated, in which 
we found a systematic decrease of the Kosterlitz-Thouless (KT) transition temperatures of both XY 
and six-state clock models as the concentration of dilution decreases. For the six-state clock model, 
a second KT transition at lower temperature was observed. The KT transition temperatures as 
well as the decay exponent rj for each concentration of dilution are estimated. It is observed that 
the quasi long range order disappears at the concentration of dilution very close to the percolation 
threshold. The decay exponent r\ of the KT transitions calculated at each concentration indicates 
that the universality class belongs to the pure XY and clock models, analogous to the expectation of 
the Harris criterion for the irrelevance of randomness in the continuous phase transition of systems 
with non-diverging specific heat. 
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I. INTRODUCTION 

Dilution is an indispensable aspect of both theoretical 
and experimental studies due to the presence of defects 
and impurities in any real material [ij, |2( ■ The pioneering 
work by Harris on the effect of dilution on the crit- 
ical behavior of systems with continuous transition has 
stimulated many studies in the field of random systems. 
Based on the Harris criterion, it is predicted that the 
dilution will be relevant (irrelevant) if the specific heat 
exponent a of the pure system is positive (negative) , and 
becomes marginal in the case a = 0, for example, in the 
two-dimensional (2D) Ising system Q. 

It is well known that the pure 2D XY model can- 
not have a true long range order at any finite temper- 
ature due to the continuous symmetry U(l) |5j,|6|. How- 
ever, the system can exist in a quasi long range order 
(QLRO), an intermediate phase which is a topological ex- 
citation formed by vortex-antivortex pairs [7j . The num- 
ber of vortex-antivortex pairs increases with the tempera- 
ture until the system experiences the Kosterlitz-Thouless 
(KT) transition. This transition corresponds to the un- 
binding of vortex-antivortex pairs, which leads the sys- 
tem to a high-temperature disordered phase. 

The existence of QLRO in the presence of dilution is 
an interesting topic. This is related to the fact that the 
QLRO is a topological order vulnerable to the local de- 
fect or perturbation. The diluted XY models may have 
the relevance to the study of superconductivity, in par- 
ticular the interaction between vortices and the spatial 
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inhomogeneity due to the impurities. However, not so 
much attention has been given to the dilution effect on 
the KT transition. 

Quite recently, two contradictory results on the 2D site 
diluted XY model were reported 0,0- By using the 
Monte Carlo (MC) simulation with the Metropolis algo- 
rithm, Leonel et al. showed that the QLRO disappears 
before the concentration of vacant sites reaches the per- 
colation threshold p c of site dilution. On the other hand, 
performing a more extensive MC study with the Wolff 
cluster algorithm 01 i Berche et al. ;9j suggested that 
the QLRO remains up to dilute concentration very close 
to the percolation threshold p c (for the site dilution on 
the square lattice p c ~ 0.593). With these two inconsis- 
tent results, it is worth carrying out a detailed study of 
the dilution problem from a different point of view. We 
here treat a bond-diluted case; it should be made clear 
whether the QLRO remains or not up to the percolation 
threshold of bond dilution. 

The effect of the g-fold symmetry-breaking fields on the 
2D XY model has been paid attention 0, [2] ■ Treating 
clock models, where only discrete angles of the XY spins 
are allowed, is essentially equivalent to probing the g-fold 
symmetry-breaking fields. The U(l) symmetry of the XY 
model is replaced by the discrete Z„ symmetry in the q- 
state clock model. It was shown [lj] that the 2D g-state 
clock model has two phase transitions of KT type at 7\ 
and T2 (Ti < T2) for q > 4. There is an intermediate 
XY-like QLRO phase between a low-temperature ordered 
phase (T < T\) and a high-temperature disordered phase 
(T > T2). The effect of dilution on the low-temperature 
ordered phase due to the discrete symmetry of the clock 
model is another interesting subject to study. 

In this paper, we study the bond-diluted XY and six- 
state clock models (q = 6) on the square lattice. We use 



2 



the Monte Carlo method with cluster flip. For the esti- 
mator of the KT transition we use the ratio of magnetic 
correlation functions with different distances 0] . In the 
next section, we describe our model and the detail of cal- 
culation method. Then, in Sec. IIIII we shall present our 
results. Sec. lIVl is devoted to the summary and conclud- 
ing remarks. 

A part of the preliminary results of the present work 
was reported at the conference, "Statistical Physics of 
Disordered Systems and its Application" , which was held 
in July 2004, at Hayama, Japan |14j . 
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II. MODEL AND SIMULATION METHOD 



The bond-diluted XY spin model is written with the 
Hamiltonian 



H 



Si ■ Si 



(i) 



where summation is performed over the whole nearest 
neighbors pairs (ij), Si being a unit planar spin vector 
occupying the i-th site of lattice system (here we deal 
with the square lattice) , and 9i the angle associated with 
the i-th spin. For the six-state clock model, the angle 
takes discrete values, 0i = 27rg/6 with q = 0, • • • , 5. The 
quenched dilution is conveyed by the coupling follow- 
ing a distribution P(Jij) = p5(Jij — J) + (1 — p)S(Jij), 
with p being the concentration of existing bonds, or we 
can say (1 — p) is the concentration of the diluted (miss- 
ing) bonds. 

We make use of the canonical sampling MC method 
with multi-cluster spin updates (Swendsen-Wang algo- 
rithm 0|). The embedded cluster idea for continuous 
spins due to Wolff ^3 i s adopted, namely by projecting 
the planar spins into a random axis so that the Kasteleyn- 
Fortuin |16( procedure on Ising spins can be performed. 
Spins with missing bonds due to dilution are not con- 
nected in forming a cluster. 

We simulated the 2D diluted XY and six-state clock 
models on the square lattice with periodic boundary con- 
ditions. We treated both models with the linear sizes of 
L = 32, 48, 64, 80, and 96. Various bond concentrations, 
from the pure case p — 1.0 down to p = 0.55, were sim- 
ulated. For each concentration of each system size we 
treated many different realizations in order to get better 
statistics; typical number of realizations is 256, except for 
p = 0.55 where more realizations were taken into account 
to compensate highly sample-dependent results. We per- 
formed 10 4 MC steps for the equilibration and 4 x 10 4 
MC steps for the measurement. We use the reweighting 
technique to obtain the thermal average at temper- 
atures different from those at which actual simulations 
were made. 
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FIG. 1: Temperature dependence of specific heat of the di- 
luted (a) XY and (b) six-state clock models for various bond 
concentrations of system sizes L = 32, 48, 64, 80, and 96. 



III. RESULTS 

A. Specific heat 

Let us start with presenting the results of the specific 
heat for the diluted XY and six-state clock models on 
the square lattice. The specific heat per spin is defined 
as follows 



C(T) 



1 



NkT 2 



l(E 2 ) {Ef 



(2) 



where k is the Boltzmann constant. The number of spins 
and the total energy are denoted by TV and E, respec- 
tively. 

The temperature dependence of specific heat for var- 
ious bond concentrations are plotted in Fig. ^ for the 
diluted XY and six-state clock models. The temperature 
is represented in units of J/k from now on. The statisti- 
cal errors are less than the order of the width of curves. 
We see that the size dependence is small, which is typical 
for the KT transition. Single finite peaks observed in the 
diluted XY model may correspond to the existence of one 
KT transition; on the other hand there are double peaks 
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FIG. 2: Ratio of magnetic correlation functions of the diluted (a) XY and (b) six-state clock models for various bond concen- 
trations of sizes L = 32, 48, 64, 80, and 96. The temperature is represented in units of J/k. 



for the diluted clock models, which signify the existence 
of two KT transitions. In both models, the positions of 
the peaks gradually shift to the lower temperature as we 
reduce the bond concentration; the peaks become rela- 
tively flat with the decrease of p. This indicates that the 
QLRO in the system gradually fades away. However, no 
abrupt change has been observed; the change with p is 
smooth and continuous. More quantitative analysis on 
the critical behavior shall be performed on the magnetic 
correlation ratio in the following. 



B. Magnetic correlation ratio 

The critical behavior of magnetic ordering can be in- 
vestigated more precisely from the evaluation of the mag- 
netic correlation function which is defined as the follow- 



ing 



g(r) = (Si ■ S i+r 



(3) 



where r is the fixed distance between spins. Precisely, 
the distance r is a vector, but we have used a simplified 
notation. 

Consider the ratio of the correlation functions with dif- 
ferent distances. At the critical point or on the critical 



line, the correlation function g{r) for an infinite system 
decays as a power of r, 



(4) 



where D is the spatial dimension and rj the decay expo- 
nent. For a finite system in the critical region, it can be 
shown that the ratio of the correlation functions with dif- 
ferent distances has a finite-size scaling (FSS) form with 
a single scaling variable, 



9(r,t,L) 
9(r',t,L) 



(5) 



if we fix two ratios, r/L and r/r', where £ is the correla- 
tion length. 

Tomita and Okabe ^| showed that this correlation 
ratio with different distances is a very good estimator for 
the analysis of the second-order phase transition as well 
as for the KT transition. The helicity modulus and 
the Binder parameter [lij are often used in the analysis of 
the KT transition, but the correlation ratio is more effi- 
cient in the sense that corrections to FSS are smaller . 
It has been successfully applied to the study of the 2D 
fully frustrated clock model 20]. 

In the present work, we consider the ratio of magnetic 
correlation functions setting r = L/2 and r' — L/A for 
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FIG. 3: The plot of T KT (L) versus l~ 2 , with I = InbL, to 
estimate KT transition temperature of the diluted (a) XY 
and (b) six-state clock models on the square lattice for L = 
32, 48, 64, 80, and 96, with bond concentration p = 0.9. For 
the clock model, both the higher (T2) and lower (Ti) KT 
transition temperatures are estimated. The data obtained 
from different values of R are shown by different marks. 
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FIG. 4: The plot of T KT (L) versus l~ 2 , with / = InfcL, to 
estimate KT transition temperature of the diluted (a) XY 
and (b) six-state clock models on the square lattice for L — 
32, 48, 64, 80, and 96, with bond concentration p = 0.6. For 
the clock model, both the higher (T2) and lower (Ti) KT 
transition temperatures are estimated. The data obtained 
from different values of R are shown by different marks. 



two distances. Thus, we evaluate the correlation ratios 
g(L/2)/g(L/4). 



1. Estimate of KT transition temperatures 

We show the temperature dependence of magnetic cor- 
relation ratio for various bond concentrations p in Fig. [2 
for the diluted XY and six-state clock models. The data 
of different sizes for each bond concentration are plot- 
ted. The statistical errors are less than the order of the 
width of curves. Although at high temperatures, curves 
of different sizes are separated, they gradually merge as 
the temperature is decreased. At lower temperatures the 
curves of different sizes are collapsed on a single curve, 
which is the behavior of the QLRO phase, in other words, 
on the critical line. The essential feature of the diluted 
systems is the same as the pure case, which indicates 
the existence of the KT transition. The approach to 



QLRO phase is described by the scaling behavior shown 
in Eq. JSJ. For the six-state clock model, the curves 
with different sizes separate again at low enough temper- 
atures. This comes from the discrete symmetry of the 
clock model, which yields the low-temperature ordered 
phase. Thus, in addition to the high-temperature KT 
transition, a second KT transition at lower temperature 
exists for the six-state clock model. We find that the 
overall behavior of the diluted clock model is essentially 
the same as the pure clock model except that the KT 
transition temperatures decrease with the dilution. We 
shall estimate the KT transition temperatures of both 
the XY and clock models using FSS. 

With the FSS analysis based on Eq. JSJ and the KT 
form of the correlation length, £ cx exp(c/\/t), where 
t = \T — Tkt|/Tkt, we can write the L dependence of 
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TABLE I: The estimates of KT transition temperatures for 
the diluted XY and six-state clock models for various bond 
concentration p. For the clock model the lower (Ti) and 
higher (T2) KT transition temperatures are estimated. 





XY 


six-state clock 




p 


Tkt 




T 2 


1.0 


0.895(3) 


0.715(3) 


0.902(3) 


0.9 


0.747(4) 


0.605(4) 


0.727(4) 


0.8 


0.575(6) 


0.489(6) 


0.574(4) 


0.7 


0.401(6) 


0.388(6) 


0.434(6) 


0.6 


0.215(12) 


0.277(4) 


0.281(4) 


0.55 


0.076(8) 


0.198(8) 


0.204(6) 



Tkt(L) as follows 



Tkt(L) — T K t 



c 2 T, 



KT 



(ln&L) 2 



(6) 



We consider the size-dependent critical temperature 
T KT (L) that gives the constant R = g(L/2)/g(L/4). In 
Fig. |3| we show the plot of Tkt(£) as a function of l~ 2 , 
with I = ln(bL), using best-fitted parameters b and c for 
the diluted XY and six-state clock models with p = 0.9. 
The system sizes arc L = 32, 48, 64, 80, and 96. For 
the value of R, 0.600, 0.650, and 0.700 are used for the 
XY model, and Tkt(L) obtained using different R are 
represented by different marks. For the estimate of T2 of 
the diluted clock model, the value of R is set to be 0.600, 
0.650, and 0.670, whereas R is 0.975, 0.985, and 0.995 for 
T\. The data using different R are collapsed on a single 
curve, which suggests that the difference of R is absorbed 
in the i?-dependence of b. The intercepts in the vertical 
axis of Fig. El give the estimate of the KT transition tem- 
peratures. The estimate of Tkt for the diluted XY model 
with p = 0.9 is 0.747(4). The number in the parentheses 
denotes the uncertainty in the last digits. In the same 
way, we estimate the two KT transition temperatures of 
the diluted six-state clock model, T2 and Ti, as 0.727(4) 
and 0.605(4), respectively. 

We also plot the data for p = 0.6 in Fig. 21 The choice 
of fixed R is different from that for p = 0.9. The esti- 
mated KT transition temperatures are given by the in- 
tercepts in the vertical axis. The KT transition temper- 
atures for p = 0.6 are lower than those for p = 0.9; for 
the six-state clock model, the QLRO phase between Ti 
and Ti becomes narrower. 

Performing the same procedure for other bond concen- 
trations, we estimate their KT transition temperatures; 
they are tabulated in Table [I] For the six-state clock 
model, the lower (Ti) and higher (T2) KT transition tem- 
peratures are estimated. In Fig. \5\ we show the phase 
diagram of the diluted systems, which is produced from 
Table U As can be seen from the phase diagram, the KT 
transition temperatures gradually decrease with dilution, 
and the QLRO phase disappears at the bond concentra- 
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FIG. 5: Phase diagram of the diluted (a) XY and (b) six-state 
clock models. As shown, there is a systematic shift of the KT 
transition temperatures as the bond concentration decreases. 
The plot suggests that the critical dilution is close to the 
percolation threshold p c = 0.5. For the clock model, due to 
the discrete symmetry, there exists a lower KT transition, 
separating the ordered phase of low temperature from the 
intermediate QLRO. 



tion which is very close to the percolation threshold p c ; 
for the bond percolation of the square lattice, p c = 0.5. 
This result is the same as that by Berche et al. Q for the 
site dilution, but different from that by Leonel et al. 0. 

For the clock model the intermediate QLRO phase is 
suppressed to a narrow range of temperature as the di- 
luted bonds increase. The lower KT transition comes 
from the discrete symmetry. Thus, the behavior near the 
percolation threshold is similar to the case of the diluted 
Ising model. The higher KT transition temperature T2 
for the clock model becomes higher than Tkt for the XY 
model with the same p near p = 0.5, which may be re- 
lated to the stabilization effect of discrete symmetry. 



2. Decay exponent 

Next, we consider the decay exponent 77 of both XY 
and six-state clock models for each bond concentration. 




FIG. 6: Double- logarithmic plot of the magnetic correlation 
function g(L/2) versus L of the diluted (a) XY and (b) six- 
state clock models, for bond concentration p = 0.9. Here the 
slope of the best-fitted line of each corresponding R gives the 
estimate of the exponent r\. 



FIG. 7: Double-logarithmic plot of the magnetic correlation 
function g(L/2) versus L of the diluted (a) XY and (b) six- 
state clock models, for bond concentration p = 0.6. Here the 
slope of the best-fitted line of each corresponding R gives the 
estimate of the exponent r\. 



We first look at the constant value of correlation ratio 
R for different sizes and find the associate correlation 
function g(L/2). We give attention to the power-law de- 
pendence of the correlation function on the system size, 
g{L/2) ~ L~' n , expressed in Eq. with D = 2. Choos- 
ing the fixed R, we have the same temperature for dif- 
ferent sizes on the critical point or on the critical line. 
Moreover, away from the critical points the same R gives 
different temperatures for different sizes, but the correc- 
tions to the power-law behavior, Eq. QJ, is the same, 
which yields the estimate of the decay exponent rj. This 
analysis of ij was used in the study of the fully frustrated 
clock model [2(j . As an example, we consider the system 
with bond concentration p = 0.9. We plot g{L/2) versus 
L for various R's in double-logarithmic scale for the di- 
luted XY and six-state clock models in Fig. El The value 
of r\ is estimated from the slope of the best-fitted line for 
each value of constant correlation ratio. Similar plots for 
p = 0.6 are given in Fig. [7| 

The independence of thus obtained rj of each bond con- 
centration for the diluted XY and six-state clock models 



are plotted in Fig. |H1 As can be seen, the exponents 
77 of bond diluted system with various bond concentra- 
tions behave similar to those of pure case, p = 1, namely 
continuously changing with the temperature in the KT 
phase. We have plotted the data down top = 0.6. All the 
data seem to be universal, and the corrections are small 
except for larger R side of the clock model with p = 0.6. 
The deviations from the pure value becomes larger for 
p = 0.55; they are not plotted here. This comes from 
the fact that it is close to the percolation threshold and 
corrections become larger. In the renormalization-group 
language, the critical behavior is affected by another fixed 
point nearby. 

Since the estimated 77 is almost constant for smaller R, 
which corresponds to higher temperature, in Fig. |SJ the 
exponent at Xkt of the XY model and T2 of the clock 
model are estimated as 

m = 0.25(1). 

This value is consistent with that for the pure case, 
1/4=0.25. For the low-temperature side (large R) of the 
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FIG. 8: Decay exponent n as function of correlation ratio R of 
the diluted (a) XY and (b) six-state clock models for various 
concentrations. 



XY model, 77 approaches 0. Meanwhile, the estimated 77 
for the six-state clock model becomes constant for large 
R. This gives the decay exponent at the lower KT tran- 
sition temperature T\ as 

771 = 0.12(1), 

which is consistent with the pure value of 1/9 = 0.111. 
The present results suggest that the exponents associated 
with the KT transitions are universal with dilution. 



IV. CONCLUDING REMARKS 

In summary, we have investigated the bond dilution of 
the XY and six-state clock models on the square lattice 



using the canonical sampling MC method with cluster 
spin updates. We have observed the KT transition of 
higher temperature separating the intermediate QLRO 
phase from the disordered phase. Due to the discrete 
symmetry of the clock model, the lower KT transition 
has been also observed. There is a systematic decrease 
in KT transition temperatures as the concentration of 
diluted bonds increases. Our estimates of the KT tran- 
sition temperatures for each concentration, both for the 
XY and six-state clock models, are listed in Tabled from 
which the corresponding phase diagram shown in Fig. [S] 
follows. As can be seen from the phase diagram, the crit- 
ical concentration of dilution is very close to the percola- 
tion threshold which is the same as the result by Berche 
et al. 9] for the site dilution, but different from that 
by Leonel et al. 8]. The bond-diluted and site-diluted 
classical spin systems show essentially the same behav- 
ior, although there may be differences for the quantum 
spin models |2ltl22l ]. Our prelimenary results for the site- 
diluted XY and clock models give the continuous phase 
transition with respect to p. Thus, our result is in favor 
of that by Berche et al. [jj. 

The phase diagram also shows that the intermediate 
QLRO phase for the clock model is suppressed to a nar- 
row range of temperature as the diluted bonds increase. 
Our estimates of decay exponents for lower concentration 
dilution suggest that the KT transition remains unaf- 
fected by dilution, which is analogous to the expectation 
of the Harris criterion that the randomness is irrelevant 
for system with non-diverging specific heat. 

Quite recently, Sasamoto and Nishimori have studied 
the phase diagram of 2D random Z q models using the du- 
ality argument 23] . The analysis of the duality argument 
for the present model will be left to a future problem. 
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